knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
library(here)
library(broom)
library(sf) # vector, spatial data 
library(tmap) # create cool thematic maps easily
library(janitor)
library(spatstat)
library(maptools)
library(raster)
library(tidyverse)
# Read in data
oil_spill_sf <- read_sf(here("data", "oil_spill_shape", "Oil_Spill_Incident_Tracking_[ds394].shp")) %>%
  select(date = DATEOFINCI, county = LOCALECOUN)

ca_sf <- read_sf(here("data", "ca_counties", "CA_Counties_TIGER2016.shp")) %>%
  clean_names() %>%
  select(county_name = name, land_area = aland)

# View CRS info
oil_spill_sf %>% st_crs()
## Coordinate Reference System:
##   User input: WGS 84 / Pseudo-Mercator 
##   wkt:
## PROJCRS["WGS 84 / Pseudo-Mercator",
##     BASEGEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["Popular Visualisation Pseudo-Mercator",
##         METHOD["Popular Visualisation Pseudo Mercator",
##             ID["EPSG",1024]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["False easting",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Web mapping and visualisation."],
##         AREA["World between 85.06°S and 85.06°N."],
##         BBOX[-85.06,-180,85.06,180]],
##     ID["EPSG",3857]]
ca_sf %>% st_crs()
## Coordinate Reference System:
##   User input: WGS 84 / Pseudo-Mercator 
##   wkt:
## PROJCRS["WGS 84 / Pseudo-Mercator",
##     BASEGEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["Popular Visualisation Pseudo-Mercator",
##         METHOD["Popular Visualisation Pseudo Mercator",
##             ID["EPSG",1024]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["False easting",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Web mapping and visualisation."],
##         AREA["World between 85.06°S and 85.06°N."],
##         BBOX[-85.06,-180,85.06,180]],
##     ID["EPSG",3857]]
# Looks like they are the same! Yay!
# Interactive exploratory plot using tmap
tmap_mode(mode = 'view')
tm_shape(ca_sf) +
  tm_borders(col = 'black') +
  tm_shape(oil_spill_sf) +
  tm_dots()

Wrangle to find oil spill incidents per county

# Join oil spill and california counties data sets/geometry
ca_oil_sf <- ca_sf %>%
  st_join(oil_spill_sf)

# Count number of oil spills per county (exclude counties that have NA values)
oil_counts_sf <- ca_oil_sf %>%
  group_by(county_name) %>%
  summarize(n_records = sum(!is.na(county_name))) 
# Choropleth Map
ggplot(data = oil_counts_sf) +
  geom_sf(aes(fill = n_records), color = 'white', size = 0.1) +
  scale_fill_gradientn(colors = c('lightgrey', 'darkorchid', 'navyblue')) +
  theme_void() +
  labs(fill = "Number of Oil Spills",
       title = "CA Oil Spills Per County 2008")

# Cool interactive choropleth map using tmap
tmap_mode(mode = 'view')
tm_shape(oil_counts_sf) +
  tm_borders(col = 'black') +
  tm_fill('n_records', palette = 'BuPu') 
oil_sp <- as(oil_spill_sf, 'Spatial') # Convert to object 'Spatial'
oil_ppp <- as(oil_sp, 'ppp') # Convert to spatial point pattern

ca_sp <- as(ca_sf, 'Spatial') # Convert to object 'Spatial'
ca_win <- as(ca_sp, 'owin') # this window will exclude marine oil spills, but there is still enough data to perform G Function analysis without these observations

# Combine as a point pattern object (points + window):
oil_full <- ppp(oil_ppp$x, oil_ppp$y, window = ca_win) 

plot(oil_full) # Illegal point (outside window) shows up as the plus sign

NEAREST NEIGHBOR (G FUNCTION)

r_vec <- seq(0, 10000, by = 100) # make vector containing values from 0 - 10,000 that will be used to calculate G(r)

gfunction <- envelope(oil_full, fun = Gest, r = r_vec, nsim = 20, nrank = 2) # Calculate the actual and theoretical G(r) values, using 20 simulations of CRS for the "theoretical" outcome (the processing time is very low so fewer simulations will run faster)
## Generating 20 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,  20.
## 
## Done.
gfunction_long <- gfunction %>%
  as.data.frame() %>%
  pivot_longer(cols = obs:hi, names_to = 'model', values_to = 'g_val')

ggplot(data = gfunction_long, aes(x = r, y = g_val, group = model)) +
  geom_line(aes(color = model))

This confirms clustering - our data has a greater proportion of events with nearest neighbor at smaller distances compared to a theoretical CSR scenario (model = theo).